COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	Barrow FFT
C00008 ENDMK
C⊗;
;Barrow FFT
;Here is the Barrow FFT benchmark which tests floating operations
;of various types, including flonum arrays. (ARRAYCALL FLONUM A I)
;accesses the I'th element of the FLONUM array A, where these arrays are
;0-based. (STORE (ARRAYCALL FLONUM A I) V) stores the value V in the
;I'th element of the FLONUM array A. 

;There was a fair amount of FLONUM GC's in the SAIL MacLisp run, which,
;when it needed to CORE up during GC, took 4.5 seconds of CPU time for the
;computation and 15 seconds for GC. Other configurations of memory required
;only 1.5 seconds for GC.

;Refer to this as FFT.
;			-rpg-

;;; *-*lisp*-* 
;;; From Rich Duda, by way of Harry Barrow -- 3/26/82 

(DEFUN FFT					        ;Fast Fourier Transform
  (AREAL AIMAG)                                         ;AREAL = real part
  (PROG							;AIMAG = imaginary part
   (AR AI PI I J K M N LE LE1 IP NV2 NM1 UR UI WR WI TR TI)
    (SETQ AR (GET AREAL 'ARRAY))			;Initialize
    (SETQ AI (GET AIMAG 'ARRAY))
    (SETQ PI 3.141592653589793)
    (SETQ N (CADR (ARRAYDIMS AR)))
    (SETQ N (1- N))
    (SETQ NV2 (// N 2))
    (SETQ NM1 (1- N))
    (SETQ M 0)						;Compute M = log(N)
    (SETQ I 1)
   L1 (COND
       ((< I N)(SETQ M (1+ M))(SETQ I (+ I I))(GO L1)))
    (COND ((NOT (EQUAL N (↑ 2 M)))
	   (PRINC "Error ... array size not a power of two.")
	   (READ)
	   (RETURN (TERPRI))))
    (SETQ J 1)						;Interchange elements
    (SETQ I 1)						;in bit-reversed order
   L3 (COND ((< I J)
	     (SETQ TR (ARRAYCALL FLONUM AR J))
	     (SETQ TI (ARRAYCALL FLONUM AI J))
	     (STORE (ARRAYCALL FLONUM AR J) (ARRAYCALL FLONUM AR
I))
	     (STORE (ARRAYCALL FLONUM AI J) (ARRAYCALL FLONUM AI
I))
	     (STORE (ARRAYCALL FLONUM AR I) TR)	
	     (STORE (ARRAYCALL FLONUM AI I) TI)))
    (SETQ K NV2)
   L6 (COND ((< K J) (SETQ J (- J K))(SETQ K (// K 2))(GO L6)))
    (SETQ J (+ J K))
    (SETQ I (1+ I))
    (COND ((< I N)(GO L3)))
    (DO L 1 (1+ L) (> L M)				;Loop thru stages
	(SETQ LE (↑ 2 L))
	(SETQ LE1 (// LE 2))
	(SETQ UR 1.0)
	(SETQ UI 0.0)
	(SETQ WR (COS (//$ PI (FLOAT LE1))))
	(SETQ WI (SIN (//$ PI (FLOAT LE1))))
	(DO J 1 (1+ J) (> J LE1)			;Loop thru butterflies
	    (DO I J (+ I LE) (> I N)			;Do a butterfly
		(SETQ IP (+ I LE1))
		(SETQ TR (-$ (*$ (ARRAYCALL FLONUM AR IP) UR)
			     (*$ (ARRAYCALL FLONUM AI IP) UI)))
		(SETQ TI (+$ (*$ (ARRAYCALL FLONUM AR IP) UI)
			     (*$ (ARRAYCALL FLONUM AI IP) UR)))
		(STORE (ARRAYCALL FLONUM AR IP)
		       (-$ (ARRAYCALL FLONUM AR I) TR))
		(STORE (ARRAYCALL FLONUM AI IP)
		       (-$ (ARRAYCALL FLONUM AI I) TI))
		(STORE (ARRAYCALL FLONUM AR I)
		       (+$ (ARRAYCALL FLONUM AR I) TR))
		(STORE (ARRAYCALL FLONUM AI I)
		       (+$ (ARRAYCALL FLONUM AI I) TI)))
	    (SETQ TR (-$ (*$ UR WR) (*$ UI WI)))
	    (SETQ TI (+$ (*$ UR WI) (*$ UI WR)))
	    (SETQ UR TR)
	    (SETQ UI TI)))
    (RETURN T)))



;;; Sets up the two arrays
(SETQ RE (ARRAY RE FLONUM 1025.))

(SETQ IM (ARRAY IM FLONUM 1025.))

;;; The timer which does 10 calls on FFT

(include "timer.lsp")
(timer timit 
       (do ((ntimes 0 (1+ ntimes)))
	   ((= ntimes 10.))
	   (fft 're 'im)))